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ABSTRACT 

We perform idealised numerical simulations of magnetic buoyancy instabilities in three 
dimensions, solving the equations of compressible magnetohydrodynamics in a model 
of the solar tachocline. In particular, we study the effects of including a highly sim- 
plified model of magnetic flux pumping in an upper layer ( "the convection zone" ) on 
magnetic buoyancy instabilities in a lower layer ("the upper parts of the radiative 
interior - including the tachocline"), to study these competing flux transport mech- 
anisms at the base of the convection zone. The results of the inclusion of this effect 
in numerical simulations of the buoyancy instability of both a preconceived magnetic 
slab and of a shear-generated magnetic layer are presented. In the former, we find that 
if we are in the regime that the downward pumping velocity is comparable with the 
Alfven speed of the magnetic layer, magnetic flux pumping is able to hold back the 
bulk of the magnetic field, with only small pockets of strong field able to rise into the 
upper layer. 

In simulations in which the magnetic layer is generated by shear, we find that the 
shear velocity is not necessarily required to exceed that of the pumping (therefore the 
kinetic energy of the shear is not required to exceed that of the overlying convection) , 
for strong localised pockets of magnetic field to be produced which can rise into the 
upper layer. This is because magnetic flux pumping acts to store the field below the 
interface, allowing it to be amplified both by the shear, and by vortical fluid motions, 
until pockets of field can achieve sufficient strength to rise into the upper layer. In 
addition, we find that the interface between the two layers is a natural location for 
the production of strong vertical gradients in the magnetic field. If these gradients 
are sufficiently strong to allow the development of magnetic buoyancy instabilities, 
strong shear is not necessarily required to drive them (c.f. previous work by Vasil & 
Brummell). We find that the addition of magnetic flux pumping appears to be able 
to assist shear-driven magnetic buoyancy in producing strong flux concentrations that 
can rise up into the convection zone from the radiative interior. 

Key words: hydrodynamics - MHD - magnetic fields - instabilities - Sun: magnetic 
fields - Sun: interior 



1 INTRODUCTION 

The Sun is observed to possess a large-scale predominantly 
toroidal field at the surface, which exhibits cyclic magnetic 
activity, as manifested by the sunspot cycle. Since the period 
of these cycles (approximately 11 yr) is much shorter than 

* E-mail: adrianjohnbarker@gmail.com 



the ohmic diffusion time for a primordial solar magnetic field 
(10 Gyr), it seems inescapable that we require the action of a 
dynamo which can generate magnetic fields to explain this 
behaviour. Since dynamos are generally assisted by differ- 
ential rotation, it is believed that the tachocline, a region 
of strong radial (and somewhat weaker latitudinal) shear at 
the base of the convection zone of the Sun, plays an impor- 
tant role in this process, by generating toroidal field from 
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poloidal fielcfl ( |Tobias fc Weiss|2007| ). The bulk of the solar 
toroidal field is likely to be stored just below the convection 
zone, both because the action of convective turbulence has 
been found to rapidly pump magnetic field into the nontur- 



bulent regions beneath (e.g. Spiegel & Weiss 1980 Tobias 
et al.|2001 i, and also because the short rise time of buoyant 



gree. Recent work has therefore begun to address the prob- 
lem of the generation of a toroidal magnetic layer through 
shear, together with the resulting magnetic buoyancy in- 
stabilities o f this layer (e.g. |Brummell et al 



magnetic flux tubes would pose problems for the storage of 



toroidal field at any depth in the convection zone (Parker 



19751 



2002 Cline 



et al.|[2003l |Vasil fc Brummeil||2008| hereafter VB08; |Vasil 



fc Brummeli|2009| hereafter VB09;|Silvers, Vasil, BrummelT 
fc Proctor||2009| |Silvers, Bushby fc Proctor||2009| hereafter 



Sunspots and other intense flux elements comprising ac- 
tive regions on the solar surface are thought to be produced 
by the instability of this large-scale toroidal field stored be- 
neath the convection zone. A prime candidate for the in- 
stability mechanism is magnetic buoyancy. This can be ex- 
plained in brief as the idea that a horizontal magnetic field 
B can support heavy gas above by virtue of the pressure 
(|B| 2 /2/x D ) that it exerts. If it is in thermal equilibrium, the 
gas in the region of the field will therefore rise, since it is 
lighter than its surroundings. A horizontal magnetic field 
that decreases with height can render the fluid top-heavy, 
which is liable to instability. 

Previous work has studied the instability of a field that 
decreases smoothly with height, in which case the unstable 
modes can be either two-dimensional "interchange modes" 
or three-dimensional "undular modes", of which the latter 
are usually dominant (see, for example, the following review 
articles: |Hughes fc Proctor|1988[|Hughes|2007| ). If the field is 
discontinuous and consists of a slab of horizontal magnetic 
field sandwiched by nonmagnetic gas in a convectively stable 
atmosphere, then the instability is of Rayleigh- Taylor type, 
and occurs for any strength of magnetic field (in the absence 
of diffusion) . In this case the instability is a two-dimensional 
interchange mode, which involves no bending of the field 



SBP09). 

Of most relevance to this work, VB08 and SBP09 used 
numerical simulations in Cartesian geometry to study the 
generation and subsequent instability of a horizontal mag- 
netic layer from an initially uniform vertical field. They 
found that strong velocity shear, which is hydrodynami- 
cally unstable to Kelvin-Helmholtz type shear instabilities, 
is required for magnetic gradients to be sufficient for mag- 
netic buoyancy instabilities to develop. Since the shear in 
the tachocline is believed to be much weaker, and hydro- 
dynamically stable, this result appeared to provide a prob- 
lem for the efficacy of this mechanism in producing buoyant 
flux structures. However, it is known that radiative diffusion 
(with diffusivity k) below the convection zone is much more 
efficient at transporting heat than ohmic diffusion (with dif- 
fusivity r)) is at transporting magnetic flux. In the regime 
that k 3> rj, it has long been known that double-diffusive 



effects could allow magnetic buoyancy instabilities ( Oilman 



1970 Acheson 1979 Schmitt fc Rosner 19831 by mitigat- 



ing the stabilising stratification through the action of ra- 
diative diffusion. It is therefore possible that a hydrody- 
namically stable tachocline shear can produce a magnetic 
layer that is unstable to a double-diffusive magnetic buoy- 
ancy instability. This was first confirmed by the simulations 



of Silvers, Vasil, Brummell & Proctor (20091, who used a 



similar setup to VB08, except in a parameter regime which 



lines (Cattaneo & Hughes 19881. However, the nonlinear favoured double-diffusive instabilities (see also Silvers et al 



evolution of the instability can generate three-dimensional 
arched structures that qualitatively resemble those observed 
at the surface. These arise through the interactions between 
vortices, which are primarily generated by Kelvin-Helmholtz 
instabilities (vorticity is also produced by baroclinicity, a lin- 
ear effect associated with the initial instability) driven by the 
rising magnetic "mushrooms" (Matthews, Hughes fe Proctor| 
1995 Wissink et al.|20"00 l. An important question concerns 
the scales of the rising magnetic structures. In the absence 
of diffusion the instability occurs at very small scales, and it 
is hard to see how such modes can lead to the large coherent 
structures seen at the surface. Larger length scales can be 
obtained either by using enhanced turbulent diffusion (and 
numerical computations are inevitably diffusive and show 
the same effect), or as a result of helicity in the magnetic 
field structures, for example through rotation. 

In reality, the predominant source of the solar toroidal 
field is likely to be the strong radial (and somewhat weaker 
latitudinal) shear in the tachocline. This effect results from 
the variation in the angular velocity in the solar interior, 
which can stretch any poloidal field to produce toroidal field. 
While any poloidal field is unlikely to be coherent in the 
tachoclinsQ it is certainly likely to be present to some de- 



20101. They showed that magnetic buoyancy instabilities of 



a shear-generated magnetic layer are possible in the double- 
diffusive regime, which is relevant for the Sun. In this paper 
we use the simplest model in which to investigate magnetic 
buoyancy instabilities, and do not study double diffusive 
effects, while recognising that a more complicated model 
should be considered in due course. 

Whatever the mechanism by which these structures are 
produced, they then rise into the solar convection zone. The 
earlier calculations on shear-induced buoyancy did not in 
general include the action of the turbulent convection in this 
region. It is likely that the convection zone plays a crucial 
role in the solar dynamo process by helping to return flux of 
an appropriate orientation to the shear region so that the dy- 
namo cycle can be completed. The effect of anisotropic and 
inhomogeneous turbulence on flux transport has been known 



for some time ( Drobyshevski & Yuferev 1974 Arter et al 



1982||Arter|1983||Oarloway fc Proctor|1983||Tao et al.|1998 l 

The principal effect is to transport horizontal magnetic flux 
down the gradient of turbulent intensity (Zcldovich 1957 



|Moffatt 19831. In the presence of a stable layer beneath, 
which is absent in convective turbulence, magnetic flux can 



1 This is the so-called £7 effect of mean-field electrodynamics 
| Moffatt|1978 l. 

this is because a coherent poloidal field which straddles the 



base of the convection zone would cause the differential rotation 
of the convection zone to be imprinted onto the radiative interior, 
which is not what we observe (e.g. [MacGrcgor & Charbo nneau| 
1999 though also see Rogers|20lT i. 
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be transported and stored in this layer 


Tobias et al.||1998 


Tobias et al.||2001| |Dorch & Nordlund| 


2001| |Ossendrijver| 


et al.|2002 


Kitchatinov & Riidiger 2008) 


Numerical Simula- 



tions of turbulent penetrative compressible convection show 
that magnetic flux is rapidly transported into the underlying 
stable layer on a convective (and not a diffusive) timescale 
by strong downflowing plumes (|Tobias et aL 1998 Tobias 



et al.||2001| porch & Nordlund 12001b, and that only the 



strongest field is able to counteract this effect and rise into 
the upper unstable layer. This phenomenon is robust and its 
physical basis is straightforward: in a stratified compressible 
convecting fluid there is a sharp distinction between rapidly 
falling and gently rising plumes ( |Weiss et al.|2004 l. 

The purpose of this paper is to carry out a pilot study of 
the effects of the turbulent pumping on the evolution of the 
buoyancy instability. The computations referred to above 
show clearly that such an interaction is meaningful, with 
flux being transported not only on the smallest scales of the 
convection but on much larger scales too, so that there is 
a 'mean' effect. Ideally we would wish to study the inter- 
action of magnetic buoyancy instabilities with the overlying 
convective flows in a direct numerical simulation, such as 
an extension of those presented in SBP09. However this is 
computationally challenging at the present time. Further- 
more, we wish to understand the basic physics of the in- 
teraction of the atmosphere and the growing instabilities, 
without getting distracted by the complexities of the ac- 
tual convection. Thus we choose to begin by constructing 
a model problem that captures the effect of the turbulence 
on the largest scales of the buoyancy modes by means of a 
mean-field ansatz. Clearly this simplification does not allow 
consideration of the largest scales of motion which are com- 
parable or larger than the important buoyancy modes. Nor 
does it properly treat the small scales of the buoyancy which 
may be comparable with or smaller than those of the con- 
vection. Nonetheless the simple ansatz used does in our view 
capture important elements of the interaction, and gives a 
useful guide to more detailed calculations in the future. Note 
that our approach will not be strictly accurate at the base 
of the convection zone, where convection exists on a wide 
range of scales (e.g. Miesch et al.|20 08l). Averaging over the 
largest convection cells would then be meaningful only for 
the largest scales of the buoyant magnetic field. However, 
these models are not designed to accurately simulate the 
conditions at the base of the convection zone of the Sun. 
Instead, they are designed to provide a simple phenomeno- 
logical picture of the magnetic pumping of large-scale fields 
resulting from small-scale convection, which can hopefully 
provide a complementary perspective to simulating the con- 
vection directly. 

In keeping with the philosophy of modelling processes 
at the simplest non trivial level, we represent the magnetic 
pumping via mean-field electrodynamics as the antisymmet- 
ric part of the a-tensor (e.g. Radlcr 1968, Moffatt 1983), i.e., 
otij = aSij —eijk'yk- This contributes an additional advective 
velocity to the induction equation for the mean magnetic 
field through the term V x (7 x B), where 7 is a turbu- 
lent pumping velocity, which arises from a gradient in the 
turbulent intensity of the flow. 

In this paper, we therefore study the effects of this 
mechanism on magnetic buoyancy instabilities by imple- 



menting a 7-pumping effect in numerical simulations of stan- 
dard magnetic buoyancy. 

We first extend previous calculations of the buoyancy 
instability of a preconceived magnetic slab (|Cattaneo fc| 



Hughes|1988l|Matthews, Hughes fc Proctor|1995| by inclucT 
ing the additional effect of magnetic pumping in an upper 
layer ("the convection zone"). This is designed to crudely 
mimic the addition of one of the most important effects of 
the convection zone on the nonlinear evolution of magnetic 
buoyancy instabilities of this magnetic slab. Though this is a 
very simple model, it is worthwhile to extend previous sim- 
ulations through the addition of only this effect, to study 
in detail its influence on the problem. Later on, we extend 
previous calculations in which this magnetic layer is gener- 
ated through shear acting on a vertical field, by including 
magnetic pumping in an upper layer. This should allow us to 
isolate the most important effect of an overlying convection 
zone on the nonlinear outcome of these instabilities. 

The structure of the paper is as follows. We first de- 
scribe the numerical setup of our simulations starting with a 
preconceived magnetic slab in fj2] In this section we also dis- 
cuss the simple model of magnetic pumping that we adopt. 
We then present the results of these simulations in Sj3] The 
numerical setup for the problem with shear is described in 
!|4] and the corresponding results in Sj5] Finally, in Sj6] we 
present a discussion of the results and a summary of those 
which we deem to be important. 



2 NUMERICAL MODEL (WITHOUT SHEAR) 

We adopt a Cartesian box with coordinates (x, y, z), which 
represents a plane section of the tachocline region of the 
Sun. The entire domain extends from z = (top) to 2 = d 
(bottom) in the vertical, and extends from x, y = to x, y — 
X Xt yd in the horizontal, where will be specified later. We 
identify xz as the poloidal plane (where — e z is considered to 
be the radial direction), and y as the toroidal (azimuthal) 
coordinate. An infinite slab of uniform toroidal magnetic 
field B = Boe y is placed in the region z £ [21,22], which is 
sandwiched by nonmagnetic gas. The initial state is static, 
with a linear temperature distribution T = To (1 + Oz/d), 
and is piecewise polytropic (with index m), being computed 
under the assumption that the total (gas plus magnetic) 
pressure is continuous at each magnetic interface. This is an 
unstable equilibrium configuration. If a small perturbation 
is imposed on the system, the density discontinuity at z = Zi 
is unstable to a Rayleigh- Taylor type instability, which we 
excite by adding small random thermal perturbations to the 
top of the magnetic layer. 

We assume the fluid is a perfect gas with constant shear 
viscosity /x, thermal conductivity k, magnetic diffusivity r] 
and specific heats c p and c v (which define the gas constant 
1Z — c p — Cv). The standard equations of mass, momentum, 
entropy and magnetic induction, together with the equation 
of state, can be non-dimensionalised by scaling lengths with 
the depth of the layer d, temperatures with the initial tem- 
perature at the upper surface To, densities with the initial 
density at the upper surface po, and magnetic fields with the 



magnitude of the initia l magnetic field Bo (e.g. Matthews, 



Proctor & Weiss 19951. We also use d/y/lZTo as the unit 



of time, which corresponds to the sound travel time, using 
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the isothermal sound speed at the top of the layer. Using 
this non-dimensionalisation, equations governing the tem- 
poral evolution of the density p, velocity u, temperature T 
and magnetic field B read 



dtp + V • (pu) = 



dt(pu) = -V (p + ^|B| 2 ^) + V • (FBB — puu + aCn'i 
+6(m + l)pe z 

d t T = -u • VT - (7 - 1)TV • u + 2^V 2 T 



+ C ^^^IMI 2 + fCoUI 2 
p v 2 



<9 t B = V x [u x B - CoCa-J] + G 
V-B = 0, 

where the current J = V x B, the equation of state is p - 
and the viscous stress tensor is 
2 



= djUi + diUj 



-8ijd k u k , 



(1) 

r) 

(2) 

(3) 

(4) 

(5) 
pT, 

(6) 



and G will be specified later. 

These equations contain seven dimensionless parame- 
ters, which, together with the initial and boundary con- 
ditions, completely determine the evolution of the system. 
These are the polytropic index m, the temperature gradient 
6, and the ratio of specific heats 7 = c p /c v , together with 
the Prandtl number 

<x = ^, (7) 

the ratio of magnetic to thermal diffusivity at the top of the 
layer 

VPoCp 



C = 

the dimensionless thermal diffusivity 
C K = 



(8) 



(9) 



(10) 



poCpdVlZTo 
and the dimensionless field strengh 

F=^^. 
1ZT p p 

The last quantity is related to the plasma /3, which is the 
ratio of gas to magnetic pressure, by F — 2//3. Note that 
the number of pressure scale heights in the domain is given 
by N p — (m + l)ln(l + 9). We always take 7 = |, as is 
appropriate for a monatomic ideal gas, and will from now 
on reuse the symbol 7 to represent magnetic flux pumping, 
which we will define in § |2.1| 

Eqs. [IJj5]are solved subject to the boundary conditions 

that 



u z = d z u x = 
T = 1 at z 



Z Uy 





0, d z T 



at z = 
at z 



(11) 
(12) 



These conditions represent impermeable, stress-free bound- 
aries at the top and bottom of the computational domain. 
All horizontal (magnetic and nonmagnetic) boundary con- 
ditions are periodic. The mass flux and mechanical energy 



flux thus vanish on the boundaries, and the imposed heat 
flux is the only flux of (non-magnetic) energy into and out 
of the system. The vertical magnetic boundary conditions 



B x = B v =0 at z = Q,l. 



(13) 



The field is therefore ensured to be vertical at the top and 
bottom boundaries. Note that any imposed horizontal field 
can leave the domain, since a gradient of these fields can 
exist at the boundaries. This choice of boundary conditions 
is somewhat artificial. However, the dynamics in the region 
not close to the vertical boundaries should only be weakly 
influenced by them. 

The numerical method used to solve the above 
system of equations is a parallel hybrid finite- 
difference/pseudospectral code, where spatial derivatives 
are calculated in Fourier space for the horizontal directions 
and fourth-order finite-differences in the vertical direction 
(upwind derivatives being used for the advection terms). 
Time integration is performed by an explicit third-order 
Adams-Bashforth method. The equations solved, and the 
numerical method used are discussed in more detail in 
Matthews, Proctor fc Weiss| ( |1995| and |Bushby fc Houghton] 
(20051, for example. The code has been thoroughly tested, 
particularly on problems of magnetoconvection and 
magnetic buoyancy. 

We simulate a box that is elongated in the direction of 



the vortex tube instability observed by [Matthews, Hughes] 
& Proctor ( 1995 1 is allowed to develop and produce three- 



dimensional structure. We use a spatial resolution of 128 x 
128 x 200, except where specified otherwise. 

2.1 Downward magnetic pumping 

We add the following term into Eq. [4] to represent the effects 
of turbulent pumping of the magnetic flux from the upper 
"convective" regions, of the form 

G = V x (7 x B) , (14) 
with 

7 = 2p [1 + tanh [(A^r 1 ( Zl - z)]] e z . (15) 

The vertical profile is designed to represent a region with 
a uniform nonzero value ("the convection zone") that 
smoothly goes to zero ( "the radiation zone" ) below a depth 
Zi ("the radiative-convective interface"). We plot an exam- 
ple of such a profile in Fig. [l] using a set of typical values for 
the various parameters. This is a purely downward pump- 
ing velocity, which should act to prevent magnetic field with 
horizontal strengths smaller than 



Bh I = B eq = 7„ 



F 



(16) 



from rising into the upper layer. The evolution of the system 
after the onset of buoyancy instabilities will therefore de- 
pend on the parameter (noting that we initially take B y — 1) 



M-, 



7m 



p(Zi 



B„ 



(17) 



which is an Alfvenic Mach number for the 7-pumping. If 
M 1 < 1, we would expect some of the rising field to be 
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Figure 1. Typical profile of magnetic pumping with Z{ 
0.5, (Ax;)" 1 = 30 and -y m = 0.1. 



able to overcome the downward pumping and rise into the 
upper layer. Alternatively, if M 7 > 1, we would expect the 
field to be too weak to overcome the downward pumping, 
and 7 should therefore act as a lid which will hold down 
the field, unless the field can be locally amplified sufficiently 
that |B fc | > B eq . 

The interpretation that we will adopt is that the Carte- 
sian box represents a large horizontal section at the base 
of the convection zone. In this case, the spatial scale of the 
convection cells is considered to be much smaller than the 
horizontal size of the box. We can then reasonably assume 
a separation of scales and define an average over the small- 
scale convection cells. Our vertical scale is such that the 
upper layer (z < Zi) represents a significant fraction of the 
lower half of the convection zone, and 7 can be considered a 
"mean-field" pumping effect. In this interpretation, the sim- 
ulated magnetic field, B, is always considered to be a "mean 
field" , with the horizontal length scales of the magnetic field 
being much larger than that of the (unresolved) convection. 
It is then reasonable to consider 7 to act with equal strength 
on all of the resolved horizontal scales of B in our simula- 
tion, to a first approximation. For numerical reasons, it is 
essential to include nonzero diffusivities of momentum, heat 
and magnetic field. This means that the resolved diffusive 
lengthscales in the simulation are, by construction, larger 
than the horizontal scales of the unresolved small-scale con- 
vection. We consider the actual diffusive lengthscales to be 
much smaller than those of the unresolved convection. If 
desired, the simulated diffusivities can be considered to rep- 
resent turbulent diffusivities resulting from the unresolved 
convection. However, we only include them for numerical 
reasons (though we will later include their effects in our dis- 
cussion, for completeness), since we do not set out to study 
the effect of turbulent diffusivities on magnetic buoyancy. 

Convective turbulence is also likely to pump the scalar 
fields of density and temperature, in addition to the mag- 
netic field. If the turbulence is incompressible and only 
varies in one direction, such scalar pumping vanishes (leav- 
ing only anisotropic diffusion) because the pumping velocity 
is divergence- free ( Moffatt|1983 1, and when it exists it is, in 
general, not simply related to the velocity of magnetic pump- 
ing ( |Cattaneo et al.||1988[ ). For compressible turbulence, it 
can be shown that an additional mean advection term arises 



(e.g. Vergassola fc Avellaneda||1997 l, which may not neces- 
sarily vanish when the intensity of the turbulence varies only 
with height, as in our case. However, for simplicity, and since 
this effect is not simply related to 7, we choose to neglect 
pumping of scalar fields throughout this paper. 



2.2 Parameters adopted 

Convective flows in the lower parts of the convection zone 
are likely to be highly subsonic, wit h mach numbers inferred 
to be in the range 10 -4 



10 2 (e.g. |Ossendrij ver 
et al.pOlO l. In the simulations of |Tobias et al. 



2003 



Jones 



(20011, the 



magnetic pumping effect resulting from turbulent compress- 
ible convection occurs on a convective timescale. We might 
therefore expect 7 to be subsonic, with velocities at most 
comparable to the convective flows, constraining 7 m <C 1. 



In the tachocline, f) > 
which constrains F < 1 



Tobias fc Hughes 20041, 



10 7 (e 

We must choose a much smaller 
/3 for these calculations than in reality to speed up the ini- 
tial instability. This is because the fastest growing Rayleigh- 
Taylor type mode has a growth rate which scales with /3 _1 
(Cattaneo & Hughes 1988J . We wish to study the nonlinear 
evolution of the instability in the presence of 7-pumping in 
the upper layer for a range of values of M 7 either side of 
unity. To do this we fix a value of F = 0.01 and vary j m . 
This fixes the growth rate and horizontal wavenumber of 
the initial instability, so that we can better isolate the con- 
sequences of varying the strength of the downward pump- 
ing. The parameter values adopted for these simulations are 
summarised in Table [l] 

At the base of the convection zone the diffusivities are 
ordered such that 1 » k » i) » f (Cough 2007 Jones et al. 
2010 1. We respect this ordering by choosing 1 3> Co ~S> cr and 
Ck *C 1, though we do take much larger diffusivities than 
are present in the Sun, as is required if we are to run fully 
resolved simulations with a sensible run time. The strat- 
ification that we adopt has approximately three pressure 
scale heights within the domain. This is designed to roughly 
correspond with a region straddling the base of the con- 



vection zone ( Christensen-Dalsgaard & Thompson 20071. 
Note, however, that the upper layer in which 7-pumping 
is present also has subadiabatic stratification (since we fix 
m > 1/(1 — 7) — 3/2), unlike in the convection zone of 
the Sun. We have also performed simulations in which the 
initial stratification changes from adiabatic to subadiabatic 
when z > 0.5, which might be more appropriate for the 
Sun. However, the results of these simulations did not dif- 
fer significantly from those with subadiabatic stratification 
throughout the box; therefore we only discuss simulations 
with a single poly tropic index in this paper, for clarity. 



3 RESULTS: INSTABILITY OF A MAGNETIC 
LAYER IN NON-MAGNETIC 
BACKGROUND WITH DOWNWARD 
PUMPING 

In this section we discuss the results of our simulations 
with a discontinuous magnetic layer initially in magneto- 
static equilibrium with its surroundings, including magnetic 
flux pumping in an upper layer. The nonlinear breakup of 
such a layer in the absence of flux pumping was studied by 
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Parameter 


Description 


Value 


rn 


Polytropic index 


1.6 


6 


Thermal stratification 


2.0 


C K 


Thermal diffusivity 


0.01 


Co 


Magnetic diffusivity 


0.01 


£7 


Prandtl number 


0.005 


F 


Magnitude of magnetic energy 


0.01 


7m 


Magnetic pumping strength 


Various 


z» 


Bottom of pumping layer 


0.5 


(Azi) -1 


Width of transition region 


30.0 


21,22 


Top & bottom of magnetic slab 


0.6, 0.8 




Box horizontal aspect ratio 


1,4 



However, localised pockets of magnetic field with strengths 
satisfying |Bh| > B eq can still be produced. Several mech- 
anisms are responsible for this. One is a nonlinear effect, 
in which rising pockets of magnetic field generate vortices, 
therefore the region in the vicinity of z% is subject to com- 
plicated interactions between them. In some cases, these in- 
teractions are able to concentrate magnetic flux below the 
interface, primarily horizontally, to produce localised pock- 
ets of strong field (we later illustrate this in a simulation with 
shear in Fig. 11 1. Two other effects which can locally concen- 



Table 1. Parameter values for the discontinuous field cases. 



Cattaneo & Hughes 


( 1988 1 in 2D, and Matthews, Hughes & 


Proctor ( 1995 \ and 


Wissink et al. (20001 in 3D. 



Since we begin with a discontinuous field, diffusion 
rapidly smears the interface of the magnetic layer, though 
this does not significantly influence the dynamics of the re- 
sulting buoyancy instabilities. The initial perturbation ki- 
netic energy decays rapidly after pushing the system from 
equilibrium, with the resulting buoyancy instability being of 
Rayleigh- Taylor type, driven by the free energy associated 
with the release of gravitational potential energy stored in 
the initial state. This instability develops by perturbing the 
upper interface of the magnetic layer, eventually resulting 
in the formation of "magnetic mushrooms". The most un- 
stable mode has a horizontal wavenumber of approximately 
four when the initial configuration is randomly perturbed. 
The local rising of the field at the top of the layer results 
in shear between field and field- free regions. This is subject 
to Kelvin-Helmholtz instabilities that produce vortices at 
these interfaces, at the sides of the mushrooms (|Cattaneo| 



& Hughes 1988). These vortices interact, and neighbouring 
vortex tubes with opposite vorticity (from adjacent mush- 
rooms) become longtitudinally unstable to an analogue of 
the (antisymmetric) Crow instability between neighbouring 
counter-rotating vortex tubes (Crow 19701. This produces 



trate the field are the vertical variation of the 7-pumping, 
which can amplify field through its nonzero divergenc^] to- 
gether with the effect of buoyancy instabilities below the 
interface driving flux upwards into the interface. A combi- 
nation of these effects can produce localised peak fields that 
satisfy |B^| > B eq , and so are able to rise, even if the layer 
as a whole is contained because Af 7 > 1. 

Rising pockets of magnetic field are ultimately sheared 
apart by interacting with the gas in the upper layer and 
do not rise far as coherent structures. Since they are pre- 
dominantly toroidal structures, they are not transversally 
supported by magnetic tension, and are therefore easily de- 
stroyed (see e.g. |Hughes et al.||1998| |Hughes fc Falle||1998 1. 
Magnetic structures that are arched in the lower layer be- 
come straightened as they rise up through the 7-transition 
region. This is because, although 7 is horizontally uniform, 
there exists a vertical transition region in which 7 increases 
with height. The peaks of the arches are therefore pushed 
downwards more strongly than the troughs, which straight- 
ens the field lines as they rise through this region. It must 
be noted that this is an artifact of our adopted profile of 
magnetic pumping, and arched magnetic structures could 
easily be produced if there were some horizontal variation 
in the strength of 7. In reality, arched magnetic structures 
could be produced either by the initial instability or its early 
nonlinear evolution, or by the action of turbulent motions in 
the convection zone on initially straight structures, with the 



latter being observed by Jouve & Brun ( 2009 1 , for example 



During this simulation, the potential energy of the ini- 
tial configuration is transformed to kinetic energy of the ini- 
tial instability. Almost immediately, the shearing motions of 
the field generate vortices through Kelvin-Helmholtz insta- 
bilities. Illustrated in Fig. [3] the integrated magnetic energy 
Volume rendering images illustrating the temporal evo- (M = Jg{^\B\ 2 )dz, where we denote a horizontal average 



three-dimensional structure in the direction of the field, in- 



ducing arching of the rising magnetic structures ( Matthews, 
Hughes fc Proctor|1995 [ 



lution of our simulation with Af 7 = 1 for |B| are presented 
111 Fig. [2] So far what we have discussed is identical to the 



evolution described in Matthews, Hughes & Proctor ( 1995 1 



and I Wissink et al. (2000), which is what we would expect 
until the field has risen far enough to reach the 7-interface. 
Once the buoyant magnetic structures reach z ~ Zi, the 
resulting evolution depends on the relative strength of the 
field compared to that which can be held down by the 7- 
pumping, i.e. whether \Bh\ > B eq , which clearly depends 
on the value of M 7 for the initial state. If M 7 <C 1 down- 
ward pumping is unable to counteract the upward transport 
of magnetic flux due to buoyancy. In this regime, the ris- 
ing magnetic structures are able to rise into the upper layer 
effectively unhindered, and on horizontal scales comparable 
to the initial instability. 

The more interesting regime is when M 7 > 1, since 
downward pumping is then efficient at holding down the 
bulk of the magnetic field, as must be the case in the Sun. 



by (X) — J J Xdxdy) initially builds until the instability 
has developed and localised strong field pockets start to rise 
into the upper layer (at t « 110). Once this occurs, the 
subsequent generation of vorticity is illustrated by the in- 
crease in kinetic energy [K — J Q (|p|u| 2 )dz) and enstrophy 
(£ = / (||V x u| 2 )<iz) within the computational domain 
for t > 110. Afterwards, the initial instability dies out and 
viscous and ohmic diffusion result in a slow decay of the 



3 This is a linear process and occurs because the effects of 7- 
pumping in the induction equation are 



9t +-yd z )B X:V 
(d t +~fdz)B z 



-B Xi yd z "/ = - (V ■ 7) B ■ e XiV 



0, 



(18) 
(19) 



for which the horizontal components contain a forcing term as 
a result of the nonzero divergence of 7. This leads to some field 
amplification in vicinity of the interface between the two layers. 
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Figure 2. Volume renderings of |B| for a simulation with M 7 = 1 at approximate times t = 46,98, 123,267 and 325, respectively. This 
illustrates the temporal evolution of the magnetic field in a simulation of the buoyancy instability of a preconceived magnetic slab with 
magnetic pumping in the upper layer. 
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Figure 3. Temporal evolution of the volume-integrated magnetic 
energy M , kinetic energy K and enstrophy E, in a simulation with 
M 1 = 1, where these quantities have been scaled as listed in the 
legend. 



Figure 5. Temporal evolution of the magnetic flux fraction con- 
tained within the lower layer the transition region <&r/<I> 
and the upper layer $(7/$, in a simulation with M 7 = 1 (solid 
lines) and M 7 = 0.2 (dashed lines). 
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Figure 4. Temporal evolution of the peak magnetic field 2 max , 
the centre of magnetic field zg , and the centre of magnetic energy 
2 H 2, in a simulation with M~, = 1. 



integrated energies. Since the influence of the (artificial) up- 
per boundary becomes more important as the simulation 
progresses, we only analyse the results until t w 400. 

It is instructive to define several measures to describe 
the spatial distribution of magnetic flux in the computa- 
tional domain, and therefore understand the relative effi- 
ciencies of magnetic buoyancy and magnetic pumping. We 
can define the vertical location of the peak magnetic field, 
the centre of magnetic field (e.g. Wissink et al.|2000 Tobias 
|et al.|200l"j ) , and the centre of magnetic energy, respectively, 
as 



z max = max {2}. 

|B| 



ZB = 



fg z(By)dZ 



(20) 



(21) 



Z B 2 



(22) 



Io(By)dz ' 

fi*(\*\ 2 )dz 

We plot these measures as a function of time in Fig.|4j where 
it can be seen that the peak field is located at a depth z « Zi. 
Note that the bulk of the field is held down in the lower layer, 
since zb > Zi- This is a result of both the interaction be- 



tween neighbouring vortices, which can act in concert to hold 
down the bulk of the field, as has been previously observed 



by Cattaneo & Hughes (19881, together with 7-pumping in 
the upper layer. Note that z B ? is located slightly higher 
than zb- This difference can occur when either B contains 
appreciable B x or B z higher up, or alternatively if the field 
contains unsigned magnetic field in this region. The latter 
can be produced by the induction of small-scale field by vor- 
tices, which is most important at the top of the transition 
region. The depth of the peak field moves around chaotically 
because magnetic field is advected by vortical fluid motions 
near the interface, though it always remains near Zj. 

We define the total magnetic flux, as well as the mag- 
netic flux contained within the lower, transition and upper 
regions as 



<E>t 



(B y )dz, 



(By)dz, 



r0.6 
/ (By)dz 
JOA 



(By)dz. 



(23) 
(24) 
(25) 
(26) 



We plot each of the last three normalised to the total <E> 
(which is a time varying quantity, primarily due to ohmic 
dissipation) at each time in Fig. [5] Magnetic pumping grad- 
ually pumps any flux out of the upper layer, competing 
against localised breakouts (and ohmic diffusion). The frac- 
tion of the magnetic flux contained in the upper layer is 
always less than approximately 5% of the total initial flux, 
which reinforces the localised nature of the breakouts. In ad- 
dition, the field is primarily of a diffuse nature in the upper 
layer because the coherence of the rising pockets of flux are 
not maintained. Note that $rj/ ( I > is slightly smaller in the 
simulation with M 1 = 0.2 than that with Af 7 = 1, which 
might seem surprising. This is because the amplification of 
field through the divergence of 7 near Zi is weaker, whereas 
the bulk of the field is held down as efficiently by vortices. 

To summarise the results of this section, we find that 
although the nonlinear evolution of buoyancy instabilities of 
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a preconceived magnetic slab can produce localised pockets 
of strong field which are able to rise up against 7-pumping, 
the bulk of the field is held down in the lower layer. This is a 
result of a combination of the interactions between vortices 
and 7-pumping in the upper layer. A combination of these 
effects might be responsible for holding down the bulk of the 
solar toroidal field, allowing only localised breakouts into the 
convection zone. 



4 NUMERICAL MODEL WITH SHEAR 

The simulations described so far have assumed a precon- 
ceived horizontal magnetic field configuration in the initial 
state. From here on, we extend our study of magnetic pump- 
ing to examine the effect of its addition on the generation 
of the horizontal magnetic layer by the action of vertical 
shear on an initially uniform vertical field, following VB08 
and SBP09. This shear is designed to mimic the radial shear 
in the tachocline. We do not consider any latitudinal (hor- 
izontal) shear because this is inferred to be much weaker 



than the radial component (e.g. Christensen-Dalsgaard & 
Thompson|2 007 ) . Our adopted vertical shear profile is 



U = 



[l+tanh[(Az s ) 1 (z 3 - z)]] 



(27) 



We choose U m < 1 so that this shear is subsonic, as is 
appropriate for the tachocline. Unfortunately, it is numeri- 
cally feasible to simulate the tachocline shear only when it is 
hydrodynamically unstable to a Kelvin-Helmholtz type in- 
stability, in which case the instabilities are standard, i.e. not 
double-diffusive, instabilities. This is because capturing the 
double-diffusive instabilities would require very high resolu- 
tion. This means that it is only possible to perform a pa- 
rameter survey when Ri = min |iV 2 / (^^) 2 | < i> which is 
unlikely to be valid in the Sun (Gough 2007} . Nevertheless, 
the resulting Kelvin-Helmholtz instabilities are suppressed 
when the induced horizontal field strength becomes suffi- 
ciently large (e.g. SBP09). This is usually well before the 
onset of any magnetic buoyancy instabilities, so it should not 
significantly influence the dynamics that we are interested in 
(except by producing inhomogeneities in the magnetic layer 
along the direction of the induced field). 

In previous work, VB08 and VB09 have studied a sim- 
ilar problem and found some difficulties in maintaining the 
initial shear profile. They used an applied stress, which is 
designed to counter viscous decay, to maintain the shear. 
However, this is not able to maintain the background shear 
in the presence of strong magnetic fields, and the shear pro- 
file that remains at the end of their simulations did not 
match that of the desired (initial) profile (see VB08 Fig. 13). 
This motivated us to consider a different approach which 
will maintain the desired shear profile and might therefore 
be better at capturing the important effects of tachocline 
shear on our problem (particularly for long duration simula- 
tions). Our approach is to decompose the velocity field into 
u = u' + Uo, and consider Uo to be steady, i.e. we consider 
the tachocline to be externally maintained. We neglect the 
back-reaction of u' on Uo, and solve Eqs.[lJ|5]for u' instead 
of the total velocity field u. This has the advantage that 



Parameter 


Description 


Value 


m 


Polytropic index 


1.6 


e 


Thermal stratification 


2.0 


C K 


Thermal diffusivity 


0.01 


Co 


Magnetic diffusivity 


0.01 


a 


Prandtl number 


0.005 


F 


Magnitude of magnetic energy 


10~ 5 


7m 


Magnetic pumping strength 


Various 


*% 


Bottom of pumping layer 


0.5 




Width of transition region 


30.0 


Um 


Magnitude of shear 


0.1 


Zs 


Location of shear 


0.75 


(Azs)- 1 


Width of shear region 


30.0 


, Ay 


Box horizontal aspect ratio 


1,4 



Table 2. Parameter values for the simulations with shear. 



the underlying shear profile is not affected by the magnetic 
buoyancy instabilities that we are interested in studying] 

We now outline the modifications to the numerical 
model outlined in [J2] which allow this problem to be sim- 
ulated. The initial state is now a (single) polytropic layer 
permeated by a uniform vertical magnetic field B = Boe z . 
In the absence of shear (Uo = 0), this is an equilibrium solu- 
tion (neglecting diffusion) . However, when our adopted shear 
profile is present, the initial state is no longer an equilibrium 
configuration because the vertical field is stretched by the 
shear to produce a horizontal field. The instability of this 
field will produce magnetic structures that buoyantly rise 
into the upper layer, where they are acted upon by down- 
ward pumping, which is implemented as in j |2.1| 

We explicitly add additional terms into Eqs. [ljj4]that 
take into account the forcing of the flow by the background 
shear Uo. These terms are 



F p = -Uo • Vp, 

F u = -UoVu'-u'-VUo, 

F T = -U • VT, 

F B = V x (Uo x B) , 



(28) 
(29) 
(30) 
(31) 



which must be added onto the right hand sides of Eqs. [I] 
[4] (in that order), where we have taken into account that 
V ■ U = 0. 



4.1 Parameters adopted 

The parameters adopted for our calculations with shear are 
outlined in Table [2] where it can be noted that most of the 
parameters remain unchanged from the model in ^2] that we 
have summarised in Table [T] 

In the presence of shear, once magnetic field gradients 
become sufficient for magnetic buoyancy instabilities to oc- 
cur, the resulting evolution will depend on the Alfvenic 
Mach number of the shear-generated magnetic layer. The 



4 |Guerrero fc Kapyla| |201l[ | used another approach in which the 
horizontally-averaged velocity profile is relaxed back to the de- 
sired shear flow over some prescribed timescale. We do not adopt 
their approach since it is likely to interfere unphysically with the 
desired dynamics of the buoyancy instabilities. 
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strength of the magnetic field in the layer will grow accord- 
ing to 

B Z U, 



B 



2Az s 



-t. 



(32) 



until such a time as the layer has spread due to the slow 
propagation of Alfven waves along the initial vertical field. 
The associated timescale for this process is t S3 Az s /(2va), 
which is the vertical AlfVen travel time across half of 
the shear region. In our units, the peak field is therefore 

max(By) S3 J7 m y^f?- This means that our parameter M 7 
from § [2] is equivalent to 



M, 



7m 



max(By) max(_Bj / ) V F 



p{Zi 



(33) 



The evolution of buoyantly unstable field which reaches the 
upper layer will therefore depend on the parameter 



M v = 



(34) 



which will play an analogous role to M 1 from § [5] (the den- 
sity ratio is approximately unity). Noting that (l/2)p7^ is 
likely to be comparable to (though smaller than) the kinetic 
energy of the convection (and remembering that 7 is not an 
actual fluid velocity in our mean-field interpretation), My is 
a measure of the ratio of the kinetic energy of the convection 
to that of the shear. If My < 1, we would expect shear to 
generate sufficiently strong horizontal magnetic fields for the 
resulting rising structures produced by magnetic buoyancy 
to be able to overcome the downward pumping. Alterna- 
tively, if Mu > 1, we would expect the generated field to 
produce buoyant structures that are too weak to overcome 
the downward pumping, unless magnetic flux can be locally 
concentrated so that |Bh > B eq . As in the simulations with- 
out shear, we study various values of My either side of unity 
by fixing U m and varying j m . 

Since the shear is designed to represent the radial shear 
in the tachocline, we choose the width of the shear re- 
gion to be much thinner than a pressure scale height, with 
Az s S3 0.03 -C 1. The magnitude of the shear is subsonic, 
with Um ~ 0.1, though this is still larger than in the Sun 
so that the evolution can be fully captured within a sensible 
run time. Note that this shear is hydrodynamically unsta- 
ble, with Ri S3 0.074 < 0.25. This means that the shear 
excites vortices (aligned in the ^-direction) through Kelvin- 
Helmholtz type instabilities. Since these instabilities are un- 
wanted, we try to partially reduce their effect by having a 
small fluid Reynolds numbei^] for the shear of Re ~ 66. In 
any case, the resulting vortices are eventually suppressed by 
the induced horizontal field. 



5 RESULTS: SHEAR PRODUCTION OF A 
MAGNETIC LAYER, ITS INSTABILITY 
AND INTERACTION WITH DOWNWARD 
PUMPING 

In this section we describe the results of a set of simulations 
with shear using the numerical model outlined in [j4] First we 

5 This is defined to be Re = u ™ Az . 




Figure 7. Horizontally integrated B y in a simulation with Mu = 
1 for multiple time intervals between t S3 15 — 150 until the field 
has built up a sufficient gradient near Zj for buoyancy instabili- 
ties to set in at this location (black lines), together with the same 
quantity after onset of buoyancy instability there, at approximate 
times t = 163 (light blue), 192 (green), 222 (red) and 266 (dark 
blue). The maximum field amplitude in the shear region is ap- 
proximately equal to Um / \/F. 



describe the temporal evolution for our fiducial case, which 
has Mu = 1. We then discuss the effects of varying Mu- 

As the simulation begins, vertical shear acting on the 
uniform vertical magnetic field generates a magnetic layer 
aligned in the (negative) y-direction (for positive U m ). In the 
initial stages this field grows in magnitude and also gradu- 
ally expands due to the vertical propagation of Alfven waves 
along the initial field. A hydrodynamic Kelvin-Helmholtz 
type instability sets in from our prescribed shear flow, which 
generates vortices aligned in the ^-direction. This occurs 
because our adopted shear is hydrodynamically unstable, 
as we explained in Q However, as the magnetic field in 
the layer grows, Lorentz forces act back on the resulting 
flow and hinder further development of the instability (e.g. 
SBP09) . It therefore eventually dies out as the field strength 
exceeds a critical value, leaving behind only some weak in- 
homogeneities in the y-direction. Volume rendering images 
illustrating the temporal evolution of our simulation with 
Mu ~ 1 for |B| are presented in Fig. [6] 

Once magnetic field gradients become sufficiently 
strong, buoyancy instabilities occur in the upper parts of 
this expanding magnetic layer. These are of two-dimensional 
interchange type (like in § [3J , which are expected to be least 
affected by the addition of an aligned shear flow (Tobias 
|fc Hug hes 2004). At the same time, vertically propagating 
Alfven waves, together with a small amount of ohmic dif- 
fusion, spread the horizontal magnetic field throughout the 
lower layer. Since the field is confined from above at the 
7-interface, it remains well confined within the lower layer 



unless its magnitude can locally exceed B eq . For this to oc- 
cur once the buoyant magnetic structures, or the bulk of the 
expanding magnetic layer, first reach the interface requires 
Mu <C 1, otherwise the bulk of the field is held down below 
the 7-interface. The maximum field strength in the shear- 
generated layer reaches B y ~ ^= ~ 35. This can be seen 
in Fig. [7] which displays the horizontally- integrated vertical 
profile of By as the field builds underneath the interface for 
t = 15 — 150 (thin black lines). The field for z > Zi builds 
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Figure 6. Volume renderings of |B| for a simulation with My = 1 at approximate times t = 15, 42, 160, 177 and 200, respectively. This 
illustrates the temporal evolution of the magnetic field in a simulation of the buoyancy instability of a shear-generated magnetic slab 
with magnetic pumping in the upper layer, using Mjj = 1. 
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up to approximately this value due to the combined action 
of shear at generating the magnetic layer and 7-pumping at 
holding down the field. 

There are several important effects resulting from the 
addition of 7-pumping in the upper layer. The simplest of 
these is to hold down the bulk of the field in the lower 
layer. If 7-pumping is not present, the bulk of the field, to- 
gether with any buoyantly unstable pockets of field, continue 
to expand throughout the upper layer. Since the Sun has 
some means of holding back field until it reaches a certain 
strength, we require some mechanism to hold down weak 
field. In our simulations magnetic flux pumping in the upper 
layer ("the convection zone") provides such a mechanism, 
with only localised pockets of strong field able to rise above 
z ~ Zi. Unlike in S|3| where the interactions between vortices 
are also capable of holding down the bulk of the magnetic 
field, in these simulations the magnetic layer would continue 
to expand because of the slow vertical propagation of Alfven 
waves along the initial field. We therefore require 7-pumping 
to hold down the bulk of the field in these simulations. 

Another important effect of 7-pumping is to produce 
strong vertical magnetic field gradients in the vicinity of Zi 
(see Fig. [7j>, which induces buoyancy instabilities at this lo- 
cation. These result from the vertical gradients of 7, and 
provide an additional way to produce strong field gradients, 
and therefore induce (in our case non-diffusive) buoyancy in- 
stabilities, without requiring strong (i.e. hydrodynamically 
unstable) shear (e.g. VB08). Once the upper interface of the 
magnetic layer is perturbed by the instability, the impor- 
tant factor is the relative strength of the field in the rising 
magnetic structures to that which can be held down by the 
7-pumping, i.e. whether \Ru\ > B eq . 

As in S|3] rising pockets of magnetic field generate vor- 
tices, so the vicinity of Zi is subject to complicated interac- 
tions between them. In some cases, the resulting fluid mo- 
tions are able to concentrate the field horizontally below the 
interface into localised pockets of strong field (an example of 
this is plotted later in Fig. 111. Note that the spatial scales 



of these magnetic structures are not necessarily the same 
as those of the initial buoyancy instability. A combination 
of the above effects, together with the continual forcing of 
the horizontal layer by the shear, can produce localised peak 
fields that satisfy |Bj,| > B eq . Such strong pockets can rise 
into the upper layer. As before, rising flux structures are 
eventually sheared apart by shear interactions with the gas 
in the upper layer, and do not rise far as coherent structures. 

As is illustrated in Fig. [8] the magnetic energy builds 
until it reaches its peak value just before the field below 
the interface has been sufficiently concentrated to be able 
to break out into the upper layer. Once buoyancy instabili- 
ties occur, the resulting local shear excites vortices, thereby 
increasing K and £ . Since we are constantly forcing the sys- 
tem through the steady shear flow Uo, the magnetic energy 
oscillates about an approximately constant value, and does 
not appreciably decay throughout the simulation. We anal- 
yse the results until t « 400, after which the influence of the 
upper boundary becomes important. 

We plot the flux fractions contained in each layer nor- 
malised to the total $ (which is a time varying quantity due 
to shear production and ohmic dissipation) at each time in 
Fig. [9] This figure shows that the magnetic flux is initially 
produced and contained within the lower layer until the field 
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Figure 8. Temporal evolution of the volume-integrated magnetic 
energy M, kinetic energy K and enstrophy S, in a simulation with 
shear and My = 1, where these quantities have been scaled as 
listed in the legend. 
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Figure 9. Temporal evolution of the magnetic flux fraction con- 
tained within the lower layer the transition region $y/<E> 
and the upper layer &u/&, in a simulation with M 7 = 1 (solid 
lines) and M 7 = 0.2 (dashed lines). 



has expanded and reached z ~ Zi. Afterwards, magnetic flux 
is primarily stored in the lower layer and transition region, 
with only a few localised outbursts transporting flux into the 
upper layer. We stress the localised nature of these pockets 
of field that rise into the upper layer, since the bulk of the 
field is well confined within the lower layer (i.e. \<&u\ •C |$| 
at all times). We also plot the various measures of the distri- 
bution of magnetic field defined by Eqs. |20(|22| as a function 
of time in Fig. |10| where it can be seen that they are each 
located within the transition region. This is the case even 
when strong localised breakouts are protruding flux into the 
upper layer. Note that again z B i is located slightly higher 
than zb- The peak magnetic field always remains at z ~ Zj. 



5.1 Variation of Mu 

We have varied Mu and looked at various values either side 
of unity to study the evolution in cases where the kinetic 
energy of the shear does and does not exceed that associ- 
ated with the downward pumping. In the Mrj > 1 regime, 
we might expect the shear-generated field to be unable to 
reach strengths sufficient for buoyant rise into the upper 
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Figure 11. Contour plot of |B|, together with velocity vectors at equally spaced points in a region of the aiz-plane, in a simulation with 
Mjj = 2 at approximate times t = 166 (left) and 169 (right). This illustrates the concentration of magnetic flux in the transition region 
by vortical fluid motions, which can produce localised pockets of strong magnetic field with strengths sufficient to rise into the upper 
layer. 
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the maximum value of Mu that we can simulate to Mu < 5, 
unless we either increase the resolution or reduce the initial 
values of j m and U m , which both require greater computa- 
tional resources. Nevertheless, simulations with 1 < Mu < 5 
indicate that it is possible in this regime for a combination 
of 7-pumping at holding down (and amplifying) the field, 
and vortical fluid motions at concentrating magnetic flux, 
to produce pockets with sufficient field strength to be able 
to rise into the upper layer. 
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Figure 10. Temporal evolution of the peak magnetic field Zmscx, 
the centre of magnetic field zg, and the centre of magnetic energy 
z B 2 for a simulation with My = 1. 



layer. However, this neglects the concentration of horizon- 
tal magnetic flux by vortical fluid motions, produced in the 
nonlinear stages of evolution of the system. We have ob- 
served the field to be locally concentrated and to rise into 
the upper layer in localised breakouts even when Mu = 4, 
i.e. the kinetic energy of the shear is significantly less than 
(l/2)p7^, and therefore also smaller than the kinetic energy 
of the convection, by more than an order of magnitude. In 
Fig. 1 1 1 1 we show an example of the concentration of field in 
the transition region by vortical fluid motions in the vicinity 
of Zi in a simulation with Mu — 2. A localised pocket of field 
with a strength Bj, > 100 ~ B eq is produced, which then 
breaks out into the upper layer. 

Note that in these simulations 7 m is not as strongly 
subsonic as it is likely to be in the tachocline. The plasma 
/3 required for a pocket of magnetic field to rise is there- 
fore much smaller than we expect in reality, and for the 
strongest field /3 <C 1. The magnetic pressure is therefore 
able to evacuate the gas within a grid cell, and this becomes 
particularly important in the strongest field pockets pro- 
duced when Mu > 1. This causes the timestep to decrease 
towards zero to satisfy the Courant-Friedrichs-Lewy stabil- 
ity constraint, and the numerical code to fail. This limits 



In this paper we have studied the interaction between mag- 
netic buoyancy instabilities and magnetic flux pumping in 
simplified numerical models of the solar tachocline. We have 
adopted an idealised "mean-field" model of magnetic flux 
pumping, in which a spatially uniform and temporally con- 
stant downward advective velocity 7 for the magnetic field 
is added into the induction equation in an upper layer, with 
no flux pumping present in the lower layer. This situation 
is designed to crudely represent the effects of magnetic flux 
pumping in the lower parts of the convection zone, which 
overly a stable radiative region containing a layer of buoy- 
antly unstable toroidal magnetic field. 

We first studied the addition of a simple 7-effect into 
simulations of the instability of a preconceived horizontal 
slab of magnetic field, exte nding the calculations of |Cat- 
Itaneo & Hughes ( 1988 1 and Matthews, Hughes & Proctor 



(1995). In this problem, the initial configuration is a mag- 
netostatic equilibrium, perturbed by small thermal pertur- 
bations, which induce Rayleigh- Taylor type instabilities at 
the top of the magnetic layer. The resulting magnetic mush- 
rooms rise until they reach the pumping layer. The effect of 
this layer on the resulting evolution depends on the ratio of 
the downward pumping velocity to the Alfven speed of the 
magnetic field, that we denote M 1 . When M 1 > 1, the mag- 
netic flux pumping effectively holds down the bulk of the 
field, only allowing localised pockets of strong field to rise, 
which have been concentrated by vortical fluid motions. 

Since the Rayleigh- Taylor type instabilities of a toroidal 
magnetic layer do not have a weak field cut off, and occur for 
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any field strength (in the absence of diffusion), there must 
exist a mechanism to hold down the field until it can exceed 



a critical strength (Hughes 20071. Magnetic flux pumping 
can provide one solution to this problem, since its addi- 
tion immediately prevents fields weaker than equipartition 
strength from rising into the convection zone. If M 7 > 1 
in reality, then magnetic flux pumping can explain why the 
bulk of the field is stored in the radiative interior, with only 
localised pockets of strong field able to rise through the con- 
vection zone towards the surface. It is interesting to note 
that our simulations in this regime show that the instabil- 
ity of a uniform initial field lying underneath a layer with 
uniform downward magnetic flux pumping can produce lo- 
calised clumps of field that rise some distance into the upper 
layer. 

We also studied the addition of a 7-effect into simula- 
tions of the instability of a shear-generated magnetic layer, 
continuing calculations in the spirit of VB08 and SBP09. 
In this problem, we consider radial tachocline shear to in- 
duce a toroidal magnetic layer, which then becomes buoy- 
antly unstable. The effect of magnetic flux pumping on this 
problem has several important contributions. One is to pro- 
duce strong magnetic field gradients near the interface of 
the pumping layer. This strongly enhances the likelihood 
of buoyancy instabilities occurring in this region, and is 
therefore one method of inducing such instabilities when the 
toroidal magnetic layer is forced by a weak, hydrodynami- 
cally stable, tachocline shear. We do not, therefore, require 
strong shear for magnetic buoyancy instabilities to be ex- 
cited (c.f. VB08; VB09). 

The evolution in cases with shear depends on the ra- 
tio of the pumping velocity to the shear velocity, that we 
denote Mu ■ One interesting result of these simulations is 
that even in the case in which Mu > 1 the shear was able 
to produce localised rising pockets of field. This is interest- 
ing because the tachocline shear is probably maintained at 
a level such that it has a similar mean kinetic energy den- 
sity as the convection. Therefore, the toroidal field that is 
produced by this shear might be expected to have at most 
equipartition strength with the convective downflows. Since 
we have observed localised pockets of magnetic flux to rise 
in the regime with Mu > 1, because field is amplified by the 
combined action of concentration by vortical fluid motions, 
shear at forcing the layer and 7-pumping at holding it down 
(and amplifying it through its non-zero divergence), this in- 
dicates that the shear is not necessarily required to be more 
energetic than the convection for superequipartition fields 
to be produced by magnetic buoyancy instabilities. Some- 
what paradoxically, magnetic flux pumping (or rather its 
strong radial gradients near the tachocline) may in fact be 
an essential ingredient in producing localised pockets of su- 
perequipartition field that are able to rise up through the 
convection zone and to the solar surface. 

One important problem, which we have not attempted 
to address here, is how does the instability produce fields 
that are sufficiently helical for the resulting magnetic struc- 
tures to survive their passage through the convection zone? 
The solution to this problem will require consideration of 
initial fields that are spatially inhomogeneous and not uni- 
directional. In addition, we have modelled the effects of mag- 
netic flux pumping in the crudest possible way. Indeed, this 
calculation does not in any sense aim to be the last word on 



the matter. Rather, it should be seen as a pilot project which 
has the limited aim of establishing the efficacy of a mecha- 
nism for flux concentration. Now that this mechanism has 
been shown to have validity the next step is a fully resolved 
calculation for more general initial fields and a turbulent 
convection zone in which the effect can be put on a proper 
quantitative footing. This work is presently in progress. 
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